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I. INTRODUCTION 

The harmonic oscillator concept occupies a central position in science and engineering due to its simplicity and exact 
solubility in both classical and quantum descriptions. Today it appears in mechanics, electromagnetism, electronics, 
optics, acoustics, astronomy, nuclear theory and so on. In the quantum approach the harmonic oscillator is one of the 
exactly solvable problems studied in detail due to its considerable physical interest and applicability. The problem of 
a one-dimensional harmonic oscillator is also one of the important paradigms leading to the enormous applications in 
a wide range of the modern physics [lj. 

There are at least three logically autonomous alternative paths to quantization. The first is the standard one utilizing 
operators in Hilbert space, developed by Heisenberg, Schrodinger, Dirac, and others in the 1920s. The second one 
relies on path integrals, and was conceived by Dirac and constructed by Feynman [2]. The third one is the phase-space 
formulation. It is based on Wigner Ys (1932) quasi-distribution function and WeylYs (1927) correspondence between 
ordinary c- number functions in phase space and quantum mechanical operators in Hilbert space |3J. This complete 
formulation is based on the Wigner function (WF), which is a quasiprobability distribution function in phase-space. 

Wigner Ys quasi-probability distribution function in phase-space is a special (Weyl- Wigner) representation of the 
density matrix. It has been useful in describing transport properties, quantum optics, nuclear physics, quantum com- 
puting, decoherence and chaos. It furnishes a third, alternative, formulation of quantum mechanics, independent of 
the conventional Hilbert space, or path integral formulations. In this logically complete and self-standing formulation, 
one need not choose sides between coordinate or momentum space. It works in full phase-space, accommodating the 
uncertainty principle; and it offers unique insights into the classical limit of quantum theory: The variables (observ- 
ables) in this formulation are c-number functions in phase space instead of operators, with the same interpretation as 
their classical counterparts, but are composed together in novel algebraic ways. 

The simplest and perhaps most straightforward generalization of nonrelativistic quantum theory towards the inclu- 
sion of relativistic kinematics leads to Hamiltonians H that involve the relativistic kinetic energy, or relativistically 
covariant form of the free energy, of a particle of mass m and momentum p, given by the square-root operator and 
a coordinate-dependent static interaction potential V(x), so H(jp,x) = \/p 2 c 2 + m 2 c 4 + V(x.) (c is velocity of light 
in free space). The eigenvalue equation of this Hamiltonian is usually called the Yspinless Salpeter equation. Y It 
may be regarded as a well-defined approximation to the BethelJSalpeter formalism [4] for the description of bound 
states within relativistic quantum field theories, obtained when assuming that all bound-state constituents interact 
instantaneously and propagate like free particles |5J. Among others, it yields semi-relativistic descriptions of hadrons 
as bound states of quarks [6j|7]. 

We study rigorously the semi-relativistic quantum ID harmonic oscillator described by the Hamiltonian operator 
composed of the relativistic kinetic energy and a static harmonic potential using Wigner formulation of quantum 
mechanics. In this research we have considered a time evolution of such system, notably the time evolution of the 
momentum, coordinate and velocity distributions, average values of their quantum operators and relativistic time 
dilation. We used both Monte-Carlo procedure and method of molecular dynamic for numerical solution of this 
problem. 
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II. WIGNER - LIOVILLE EQUATION 

Integral form of the Wigner - Lioville equation. The most conventional formulation of quantum mechanics 
is description of system's dynamics with complex wave function. This function of particle with mass m and moving 
in potential field V(x) (we consider one dimension) satisfies the Schroedinger equation: 

^H ^ ^=Hn q ,t) (1) 

, where hamiltonian is H = p 2 /2m + V(q) and the initial condition is *&(q, 0) = ^o(q) 

More general description of quantum systems is given in terms of the density matrix p(q,q f ;t) (in coordinate 
representation), which has the following form in case of pure state: p(q,q'',t) = ^(g, t)ty*(q', t). Evolution equation 
for density matrix is 

^=(5-^)p W ;t), (2) 

where H acts on coordinate while H' acts on coordinate q' . The initial condition for this equation has form 

P(q,q';0) = Po(q,q') = Vo(q)n(q')- 

In the Wigner representation of quantum mechanics we use joint distribution of quasiprobability for momentum 
and coordinate; it is called Wigner function |3J. Wigner function is defined as Fourier transform of density matrix on 
difference variable £ = q — q' : while center varible is x = {q + q')/2\ 

oo 

W(x, P5 *) = ^ / P(* - x + Z/2)j*l h d£. (3) 

-co 

Evolution of quantum system in Wigner representation is describing by Wigner - Liouwille equation |3J: 

dW(x, p\ t) dH(p, x) dW{x, p- 1) dH(p, x) dW{x, p; t) 



dt dp dx dx dp 

1 



ih 



ih d \ ( ih d\ N d 

v ^ + ^)- v { x -^d P )- lhv ^ ) d P 



W(x,p;t), (4) 



with initial condition W(x,p;0) = Wq(x,p) defined by ([3| at t = 0. When the related to Hamiltonian the classical 
Hamilton's function has the form H(p, q) = p 2 /2m + V(x), the partial derivatives are: dH/dx = V^(x), dH/dp = p/m. 

In quantum mechanics Wigner functions are real valued but altering sign analog of the probabilistic joint p and x 
distributions in classical mechanics. This is supported by its general properties |3J: a) density of probability in mo- 

co oo 

mentum space is W(p) = j W(x,p)dx; b) density of probability in configuration space is W(x) = j W(x,p)dp . 

— CO —CO 

In additition W(x,p) is bilinear in wave function *$f(x,t) or in \£(p, t) in momentum representation; 
One can rewrite evolution equation Q in the integral form [3j [8] : 



W(x,p;t) = J G(p,x,t;p ,XQ,0)Wo(xQ,p )dpQdxQ + 

dr dp r dx r G (p,x,t]p T ,x T ,r) / ds W (x r ,p r — s, r) u (s, x T ) (5) 

JO J J -oo 

with lj(sjx) = F (x) ^jft _j_ _4_ J dqV (x — q) sin (2sq). Here F (x) — —dV(x)/dx is the classical force. Here 
G(p, x, t;p T , x T , r) is the Green's function for classical Liouwille equation 

G(p,X,t]p T ,X T ,r) = S(p - p(t]T,p T ,X T ))6(x - x(t]T,p T ,X r )), (6) 

where x are solutions of Hamilton's equations 

dp dH(p,x) dx dH(p,x) 

~dt = dx ' dt = dp ' ^ 
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with initial conditions p(r; r,p r , x r ) = p T , x(r; r,p r , x T ) = ab- 
solution of equation ([5| can be written in the form of the iterative series, which have the following interpretation. 
The first term of the series (the first term in the r.h.s of the Eqs |5| is equal to the sum of the contributions of the 
virtual classical trajectories defined by the dynamical Hamilton's equations ([7]). The contribution of the each virtual 
trajectory is equal to the value of the initial Wigner function Wo(xq^pq) taken at initial point xq,pq. Next terms of 
the iterative series are equal to the sum of the contributions of virtual trajectories consisting of segments of classical 
trajectories, separated by 'jumps' in momentum. The term's number in iterative series is equal to number of the 
momentum 'jumps' as it follows from convolution structure of integral term in equation (J5|. In the classical limit 
(ft — ?► 0) the force term in uj cancels the last term and only the first term of iterative series gives the main contribution 
to solution of the Wigner-Liouville equation. In the classical limit Eqs. Q, Q are reduced to the classical Liouville 
equations. 

Restrictions on initial condition. To find solutions of the integral equation ([5| the initial function Wo(x,p) 
have to be taken according to the definition (J3|. However this definition imposes the certain restriction on the choice 
of possible functions in phase space. This can be easily illustrated for harmonic oscillator (particle in potential field 
V(x) = Vo + V\x + V2X 2 ). In this case the Wigner-Liouwille equation Q is reduced to the form 

dW(x,p- 1) p dW{x,p-t) _ ^ dW(x,p- 1) _ 
at m ox op 

and coincides with the classical Liouwille equation. However solution of quantum harmonic oscillator is totally 
distinguished from classical one. Consequently, quantum and classical solutions of this task distinguish from each 
other in choice of the initial condition Wo(x,p). So we need additional condition to choose classical or quantum 
solution. For density matrix of pure state this condition can be formulated as follows 

CO 

d 2 . . d 2 [ TTr/ xi+X2 



8x18x2 



In p( Xl , x 2 ; 0) =0, or in Wigner representation [3J In / W ( Xl ^ X2 ,p)e ip(xi - X2)/h dp. (9) 



In |9J one can find another condition of the choice Wo(x,p), which is equivalent to this one. One of the physical 
interpretation of this additional conditions is connected with requirement for mean-square coordinate and momentum 
deviations to satisfy Heisenberg's principle of uncertainty at time equal to zero. 

Physical quantities. To calculate the average value of physical quantity < A > corresponded to quantum operator 
A the Weil's symbol A(p,x) has to be introduced by expression |3J: 

CO 

A(p,x)= J expip£/h(x - | A x+ d£ then < A >= J ' W(x,p)A(p,x) dpdx, (10) 

-co 

where W(x,p) is Wigner function. 

Non-relativistic harmonic oscillator. Non-relativistic harmonic oscillator with mass m and circular frequency 
uj has Hamiltonian 

A = £ + =^; (n, 

Wigner-Liouville equation Q in case of harmonic potential V(x) = muj 2 x 2 /2 has a simple form: 

dW(x,p;t) p dW(x,p;t) 2 dW(x,p;t) 

^ 1 ^ muj x =0, (12) 

ot m ox op 

As the initial condition for equations (12) we will consider a coherent state of harmonic oscillator given by the wave 
function [10J : tf(a;,0) = (ff) 1/4 exp |% - V^pf j exp | _ im j . Here p,x are average values of momentum and 
coordinate at the initial moment t = 0. According to definition (J3|, the initial Wigner function is defined by: 



muj(x — x) 2 (p — p)' 



W ( X , P ) = - exp I ^ - ^ }. (13) 



Solution of equation ( (12] ) as well as the related integral equation ([5| can be written in the form 
W(x,p\t) = J (S(p - p(t]0,po,xo))S(x - x(t]0,po,xo)))Wo(xo,po) dpodxo = 

■}• 



1 [ m uj r , _p . o 1 ^ o I 

exp < —[x — [x cos out H sinujt)\ — [p — [pcosuot — mujxsinuut)\ >. (14) 



yj2n(K/2mui) \ ft ' rnuj hmuj 
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where p(t), x(t) are virtual trajectories defined by the Hamilton's function: H(p, x) = ^ + mu; ^ x2 and the Hamilton's 
equations (FtI) with initial conditions p(0) =Po, x(0) = xq: p(t) = po cos ojt — mojx o sin ujt^x(t) = sinuot + #o cosu;t. 
Wit! 



From (14) it follows that in a coherent state of harmonic oscillator the average momentum and coordinate satisfy 



classical law of motion < p(t) >= p cos ut— mux sinuot, < x(t) >= xcosut+^sinut and the standard deviations 
of momentum and coordinate are constant: < Sp 2 >= < Sx 2 >= While Heisenberg's formula of 

uncertainty has it's minimum: < Sp 2 >< Sx 2 >= h 2 /4. Average energy is constant equal to E = ^f- + + mu; 2 x . 
More general case of the composite states slightly distinguished from pure coherent state is considered in [8]. 

Semi-relativistic harmonic oscillator. The Hamiltonian of the semi-relativistic harmonic oscillator has the 
form: 

9^9 

mu x 



H = ^p 2 c 2 + m 2 c A + (15) 



Wigner function has to be a solution of the Wigner - Lioville equation Q jllj : 

dW(x,p;t) pc 2 dW(x,p;t) 2 dW(x,p;t) 



mu x — — — = 0, (16) 



dt y/p^c 2 + ra 2 c 4 <9x "' <%> 

with initial condition W(x,p;0) = Wo(x,p) and the Hamilton's function: H(p,x) = ^/p 2 c 2 + ra 2 c' 
As before solution of this equation looks like: 



,4 j_ mo; x 



W(x,p\t) = j G(p,x,t;po,xo,0)Wo(x ,po)dpodx , (17) 

where in the Green's function the virtual trajectories p(t; to,po, #o) and x(t; £o?Po> #o) are solutions of Hamilton's 
equations 0. When c — > oo the Hamilton's function is equal to its non-relativistic limit (with rest energy term mc 2 ) 



Period of oscillations of the virtual trajectories versus energy have been found in [12J. Indeed from ^ one can 
obtained an equation for p(t) and it's first integral: 



d 2 p , 2 P n ^ d P X 2 



■ muj 2 c 



dt 2 ^2 + m 2 c 2 5 



= 0, Ci = ^ J + 2mu 2 ^p 2 c 2 +m 2 c±, (18) 



where Ci = 2mu 2 E, energy E = \/pqC 2 + m 2 c 4 + mu 2 x\j2 has to be founded from initial conditions po an d #o- So, 
after integration of this equation, the period of oscillations depends on energy of the trajectory 



^y^ + 1 ^(E-mc 2 \ 4V2 T (E-mc 2 



TW = ^ B VE + m?/ ZTSrT\ K IcTm?' ' (19) 

Y rnc 2 * 

Here i^(^) and E(z) are elliptical integrals of the first and second kind [13J: 



K(z)= - " , yJl-zsm 2 ^ (20) 

Jo v 1 — z sin Jo 

This dependence is represented by the right panel of the Fig. [5] 

Numerical simulation. We are going to obtain and compare results for quantum non-relativistic and semi- 
relativistic harmonic oscillators with the same initial Wigner function, which corresponds to a coherent state of 
non-relativistic harmonic oscillator (for simplicity we chose the average momentum po = see (p~3|) ) : 

W (x ,Po) = \exp 



mu(x - xp) 2 pi 
h hrnuj 



As we has mentioned before, an evolution of semi-relativistic oscillator is described by the expression (17). To 



consider the time evolution of oscillator the following numerical procedure has been used. The initial Wigner function 



(21) was considered as probabilistic distribution of points (po?#o) i n phase space. To sample these points we used 



Monte Carlo procedure. Each (po?#o) point was considered as the initial point of the virtual dynamic trajectory, 
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described by Hamilton's equations To solve these equations we used molecular dynamics method. Distribution 
of virtual trajectories in phase space allow to obtain Wigner function W(x,p]i) at any time t. 

Average values of general quantum operators can be obtained by calculation of the time dependences of the Weyl's 
symbol of operators along the virtual trajectories and averaging over ensemble of all trajectories. For average values 
of momentum, coordinate, energy and mean -square values of the momentum and coordinate these averaging have 
been done at each time in considered time evolution interval. 

In our simulations we have generated 10 5 of virtual trajectories. For simulation time dynamics we used implicit 
finite- difference scheme with centering [14J : 



Pm(i+1) Pm(i) 



c m(i+l) 



x m(i+l) x m(i) + 



Pm(i+1) 



Pm(i) 



xj z(z + p m{i+1) 2 ) yj 'z(z + Pm(i) 2 ) 



(22) 



Here r = 0.01 is a time step; this system of algebraic equations was mainly solved by method of simple iterations. 



For numerical calculations we used the following system of unities: p = ^-\fzp. 



x = 



t lb tm, where p n 



x m: t m are values in scheme (22) for h = 1, c = 1 However, all results below are written in usual physical units. 

Physical dimensionless parameter z = mc 2 /ftw defines "degree of relativism" of the oscillator, when z » 1 
relativist ic effects almost disappear. 



III. NUMERICAL RESULTS 



Time evolution of the distribution functions. Let us consider four different oscillators with equal parameter 
lj, but with different masses m, so the related parameters z = mc 2 /hu are equal to 100, 10, 1, 0.1 (lines 1, 2, 3 and 4 
respectively on figures below). The time evolution of momentum, coordinate and velocity distributions are presented 
by the FigsjTJ [2j [3j The almost non relativistic oscillator relates to large value of z = 100. As it follows from Fig [I] for 
z = 100, the Gaussian shape of momentum distribution is conserved during the time evolution. The same is valid for 
coordinate distribution (not shown). Just on the contrary the shape of the momentum and coordinate distributions 
for relativistic quantum oscillator are considerably changing with time (Figsfl] [2|. Firstly, one can see significant 
distribution spreading; secondly, tails of the distributions are "drawn forward" due to the difference in the periods 
of oscillation of the virtual trajectories with different initial energies especially for larger values of momentum and 
coordinate (see (19)). Oscillations of the trajectories with larger values of energy retard by phase from oscillations 
with lower energies. This results in appearance of unexpected local maximums ('protuberances') (Figjlj top left 
panel). Let us stress that the initial W(p) and W{x) (t = 0) are the normal Gauss distributions. 

Velocity distribution W(v) has an interesting shape presented by the Fig. [3] for oscillators with z = 1 (top) and 
z = 0.1 (bottom). Due to the complicated transformation v = , pc positions of the maximum of the velocity 

yj p 2 +m 2 c 2 

distributions W(v) at initial time t = does not coincide with the position of the maximum of momentum distribution 
at p = 0. Moreover at time evolution the velocity distribution is restricted by the light speed velocity, while the 
momentum distributions has not any limits. Asymmetry of velocity distribution at the initial moment is the result of 
errors, which is related to sampling of the exponentially rare events in the tails of the momentum Gauss distributions 
function. Take notice that on Fig. [T]- Fig. [2] distributions are not normalized on unity. Its presented values show the 
number of the virtual trajectories have been counted in the vicinity of the each point on the plane related to presented 
distributions. 

Average values of quantum operators. Now we are going to consider behavior of average values of quantum 
operators in time for discussed above initial Wigner function. On the Fig. [4] one can see time dependence of average 
momentum < p(t) > and the average coordinate < x(t) >. For z = 100 (weak relativism) < p(t) > and < x(t) > 
are sinusoidal functions with period of oscillation equal to 2tt/uj. Here the average momentum and coordinate behave 
almost classically like sinusoidal trajectory corresponding to initial data p(0) = 0, x(0) = xo [ID]. To analyze the 
increasing influence of relativistic effects let us consider lines 2,3,4. Firstly, with decreasing parameter z period of 
oscillations is increasing. Secondly, oscillations are damped as the Wigner functions are spreading in the phase space. 

Fig. [5] presents the time evolution of the average velocity < v(t) >. Behavior of the average velocity is similar 
to behavior of the average momentum and coordinate. Moreover, the phases of momentum and velocity oscillations 
(Fig.[4| coincide with each other. Let us stress that due to the influence of relativistic effects (lines 2,3,4) < p > / < 
v > is not equal to mass of oscillator m. 

The products of momentum and coordinate dispersions < Sp 2 >< 5x 2 > versus time are presented on the left panel 
of Fig. 6j The Heisenberg's uncertainty principle is satisfied as < 5p 2 >< 5x 2 >> ft 2 /4,. More over these products 
have minimum at t = due to the proper choice of initial state. In case of weak relativism z = 100 (line 1) the 
product of dispersions is almost constant as the initial Wigner function is not spreading. This results to the horizontal 
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FIG. 1: Time evolution of the non normalized momentum distributions for oscillators at z — 100, z — 1 and z — 0.1. The left 
top panel - initial at t = (solid line) and at uot — 86.7 distributions (dot line) for z — 0.1 . The right top panel - the contour 
plot for z — 100, the left bottom panel - the contour plot for z = 1, the right bottom panel - the contour plot for z = 0.1. 
Small oscillations on the initial distribution at t = characterize the error of the generated Gauss momentum distribution (21 ). 





FIG. 2: The time contour plot of the coordinate non normalized distributions. Left panel - z = 1, right panel - z — 0.1. 



line presenting the < Sp 2 >< Sx 2 >. On the contrary the products of dispersions increase in time very fast for small 
z. The reason for this is that the initial Wigner function does not correspond to the coherent or eigen states of the 
semi-relativistic harmonic oscillator. 

The right panel of Fig. [6] show the dependence of average energy < E(t) > on time. As the considered oscillators 
are conservative the energies are constants. For weak relativism (lines 1, 2) average energy < E > is circa the rest 
energy mc 2 (100hw and 10hw ). Due to the strong relativistic effects (line 4) the average energy < E > is considerably 
greater then rest energy mc 2 (mc 2 = O.IHuj for z = 0.1). 

Virtual trajectories. To understand physical reasons of the different properties of non-relativistic and relativistic 
harmonic oscillators we have to consider behavior of individual virtual trajectories related to initial delta distribution 
with fixed initial momentum and coordinate 



Wo(x ,po;0) = 6(x-x )6(p-p ), 



(23) 



where (po,xq) are initial condition of the virtual trajectories. Here due to the limitation (|9| distribution Wo can not 
be considered as physical Wigner function. However this choice of the Wigner function can help in understanding 
peculiarities in the time behavior of the discussed above average values. The given choice of the Wigner function 
(23) reduces the averaging over ensemble of the virtual trajectories to the consideration of the contribution only one 
virtual trajectory with fixed initial xq,pq = like in classical case. 
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FIG. 3: Time evolution of the non normalized velocity distributions for oscillators at z — 1 (top) and z = 0.1 (bottom). The 
left top panel - initial at t = (solid line) and at uot — 86.7 distributions (dot line). The right top panel - the time contour 
plot, the left bottom panel - initial at t = (solid line) and at uot = 86.7 distributions (dot line), the right bottom panel - the 
time contour plot 



Fig. [7] presents the virtual trajectories as function of time. For z = 100 non-relativistic trajectories look like 
sinusoid curves with equal oscillation periods (trajectories 1, 2 on the left panel). Let us note that energy of the 
virtual trajectory is proportional to the square of xq and in the relativistic case (central and right panels z = 0.1) 
oscillation period strongly increases with growth of the energy of trajectory (see (fl9|)). Period of oscillations of x(ujt) 
for curve 4 is longer than for one 3 due to the larger value of energy. Moreover, curves 3 and 4 tend to become the 
zig - zag lines. 

Velocity of the trajectory 4 (right panel of the Fig. [7]) looks like rectangular wave for high energy. This figure 
illustrates that mostly the virtual trajectory has velocity v approaching the velocity of light c and only near the 
turning points v differs from this limit. Velocity of trajectory 3 for the same z but lower energy differs considerably 
from rectangular wave. 

Left and cetral panels of the Fig. [8] show non-relativistic momentum virtual trajectory versus time for z = 100 (lines 
1,2) and relativistic trajectories for z = 0.1 (lines 3,4). The time dependences of p(ujt) look practically as sinusoid 
curves. In the left panel period of momentum oscillations for trajectories 1 and 2 is almost the same. In the central 
panel oscillation period is longer for the relativistic trajectory 4 with higher energy. 

Comparison of analytical and numerical calculations of the oscillation period is presented by the right panel of the 
Fig. [8j Analytical dependence of oscillation period on energy according to the formula ( 19 ) is plotted by solid curve 



for z = 0.1. Results of numerical calculations presented by triangles agree very well with analytical dependence. 

Time dilation. One of the most important relativistic effects is the time dilation: proper time of relativistic 
particle is slower than time in the non -relativistic lab frame of reference (inertial). Now we will consider this effect for 
semi-relativistic harmonic oscillator. We begin our consideration for one virtual trajectory p(t),x(t). The relationship 
between proper time and lab time is well known for the case, when particle moves with constant velocity: 



f-f = (t-t )^l-^, (24) 

where x(t) is velocity of the particle, t — to is time interval between two events in lab frame of reference, t' — t' Q is time 
interval between this two events in the rest frame of particle. In case of oscillator, however, particle's velocity is not 



a constant and formula (24) is right only for infinitesimal time intervals: 



dt' = dt\ll - while for finite time interval t' - ^ J ^ ' 1 X ^ 



*o=/^/l-^A (25) 



to 
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FIG. 4: Dependence of average momentum < p(t) > and average coordinate < x(t) > on time. Lines 1 correspond to 
oscillator with z — 100, line 2 - z — 10, line 3 - z = 1 and line 4 - 2? = 0.1. With decreasing of z period of oscillations rises, 
while oscillations are damped. 




FIG. 5: Dependence of average velocity < v(t) > on time. Line 1 corresponds to oscillator with z = 100, line 2 - z — 10, line 
3 - z = 1 and line 4 - z = 0.1. 



Dependence on time for virtual trajectories with different energies is represented in logarithmic scale by the left panel 
of Fig. [9] for z = 0.1. Almost constant parts of the curves ((4), (5)) related to the trajectories with high energy and 
correspond to the motion with velocity of order of the speed of light c. 



In case of semi-relativistic oscillator with Hamiltonian (|15| one has to average integral in (25) over all virtual tra- 
jectories. So for quantum oscillator described by the Wigner function W(x,p;t) the time dilation for semi-relativistic 
oscillator can be calculated by formula: 



f f 1 / x(t) 2 

t' Q = dt W(x,p;t) dpdx related to the 'time dilation operator'- di! = dty 1 5—. (26) 

J J yl + (p/mc) 2 V c 

to 



In numerical calculation the averaging has been done over all virtual trajectories for each time moment t. Results are 
represented by the right panel of the Fig. [9]for oscillators with z = 100, 10,1 and 0.1. 
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FIG. 6: Left panel - dependences on time of the products of the momentum and coordinate dispersions. Right panel - 
dependence of average energy < E > on time for semi-relativistic oscillators. 
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FIG. 7: Virtual trajectories for p = 0. Left panel - non-relativistic trajectories 1 (with xq = 0.05c/o;) and 2 (with xq = 0.2c/ 'uj) 
for z = 100. Central and right panels - relativistic trajectories 3 (with xq = 1.58c/o;) and 4 (with xq = 6.3c/ uj) z = 0.1. 



IV. CONCLUSION 



In this paper we are doing the exact simulation of time evolution of semi-relativistic quantum ID harmonic oscillator. 
To solve the Wigner - Liouwille equation for such system we combine Monte-Carlo procedure and molecular dynamics 
methods. As initial Wigner quasi distribution function we have used a coherent state of appropriate non relativistic 
harmonic oscillator. We have studied the time evolution of the momentum, velocity and coordinate Wigner distribu- 
tions and average values of quantum operators. Obtained results demonstrates, that relativistic treatment results in 
the appearance of the new physical effects as opposed to non-relativistic case. Interesting is the complete changing 
of the shape of the momentum, velocity and coordinate distribution functions as well as formation of "unexpected 
protuberances". We have also calculated relativistic time dilation for oscillator, as it can be useful for consideration 
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FIG. 8: Virtual trajectories for po — 0. Left panel - non relativistic trajectories 1, 2, central panel - relativistic ones 3 and 
4. Oscillation period of relativistic trajectory 4 is longer than 3 due to the larger value of virtual energy. Right panel - period 
of oscillations of the virtual trajectories versus energy for z = 0.1. Solid line - formula (19), triangles - results of numerical 
simulations, dashed line - non-relativistic case. 
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FIG. 9: Time dilation for to = 0. Left panel - the time dilation for virtual trajectories with z — 0.1. Energies E tr of the 
virtual trajectories: 1 - O.lftw (resting particle); 2 - 0.13125/kj; 3 - 0.225/kj; 4 - 0.6/kj; 5 - 2.1/kj. Right panel - the time dilation 
for quantum semi-relativist ic oscillators. 



of the life -time of the particle bound states in the traps. 
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